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Centered finite volume methods are considered in the context of Numerical Relativity. A specific 
formulation is presented, in which third-order space accuracy is reached by using a piecewise-linear 
reconstruction. This formulation can be interpreted as an 'adaptive viscosity' modification of cen- 
tered finite difference algorithms. These points are fully confirmed by ID black-hole simulations. In 
the 3D case, evidence is found that the use of a conformal decomposition is a key ingredient for the 
robustness of black hole numerical codes. 
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I. INTRODUCTION 

In recent times, many numerical relativity groups have 
performed long-term binary-black-hole (BBH) simula- 
tions. This was a long sought goal, with the main objec- 
tive of computing gravitational wave patterns that can 
be used as templates for detection. The BBH case was 
specially relevant in this respect because it is assumed to 
be the most likely candidate for detection by the current 
ground-based interferometric facilities. This can explain 
why the focus in these simulations has been placed in 
the accurate modelling of the wave zone: the numerical 
boundaries are placed safely far away, which implies the 
use of large computational domains. Also, the ability to 
extract the gravitational wave signal from the wave zone 
evolution requires the simulation to last for quite a long 
time. These facts, together with the use of some form 
of mesh refinement in order to ensure the required ac- 
curacy, make BBH simulations very demanding for the 
computational point of view, requiring a big computing 
infrastructure. 

Black hole simulations, however, deserve some interest 
by themselves, independently of the quest for gravita- 
tional waves. One can focus for instance on the strong 
field region, which can be modelled by using modest-size 
computational domains. In this case, one must refrain 
from excising the black hole interior, although many in- 
teresting results have been obtained by using excision [l[ , 
even in cases with some matter content [2j. The conse- 
quences of this choice are well known (see Ref. 0] for a 
very clear recent example): 

• A singularity- avoidant gauge condition must be 
used in order to prevent a singularity to form in- 
side the computational domain in a finite amount 
of coordinate time. 

• This makes the lapse to collapse in the black hole 
interior zones, while keeping its initial profile in the 
outer region. 

• As a consequence, steep lapse gradients appear near 
the apparent horizon, which challenge the stability 
of the numerical algorithm. 

Most of the current BH simulations are performed with 
finite difference algorithms. Regarding space accuracy, 



the most common approach is to use a centered fourth- 
order accurate method, combined with some artificial dis- 
sipation term (Kreiss-Oliger dissipation) Q. The lead- 
ing error in the solution is precisely the artificial dissi- 
pation one, usually of fourth order. One can interpret 
this combination just as a particular third-order scheme 
with some built-in dissipation, which can be tuned by a 
single parameter. This may be a difficulty in some cases, 
where dealing with the black hole interior would require 
an amount of dissipation which can be instead too big 
for the exterior region (see for instance Ref. [H). Our 
point is that centered Finite Volume methods can pro- 
vide alternative third-order accurate algorithms in which 
the built-in dissipation is automatically adapted to the 
requirements of either the interior or exterior black hole 
regions. 

Finite Volume (FV) methods have a reputation of be- 
ing computationally expensive, a price that is not worth 
to pay for spacetime simulations, where the dynamical 
fields usually have smooth profiles. From this point of 
view, centered FV methods can provide some improve- 
ment, because the they do not require the full character- 
istic decomposition of the set of dynamical fields: only 
the values of the propagation speeds are needed [j| . 

This point can be illustrated by comparing the clas- 
sical FV techniques implemented in a previous work Q 
(hereafter referred as paper I) with the new FV meth- 
ods presented here. In paper I, the general relativistic 
analogous of the Riemann problem must be solved at ev- 
ery single interface. This implies transforming back and 
forth between the primitive variables (the ones in which 
the equations are expressed) and the characteristic ones 
(the eigenvectors of the characteristic matrix along the 
given axis). In the present paper, a simple flux formula is 
applied directly on the primitive variables, so that switch- 
ing to the characteristic ones is no longer required. The 
flux formula requires just the knowledge of the charac- 
teristic speeds, not the full decomposition. 

Another important difference is that in paper I, the 
primitive quantities where reconstructed from their aver- 
age values in a piecewise linear way, using a unique slope 
at every computational cell. Only (piecewise) second or- 
der accuracy can be achieved in this way, so that going to 
(piecewise) third order would require the use of 'piecewise 
parabolic methods' (PPM), with the corresponding com- 
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putational overload. In this paper instead we split every 
flux into two components before the piecewise-linear re- 
construction (flux-splitting approach Q). This allows us- 
ing a different slope for every flux component: this extra 
degree of freedom allows us to get (piecewise) third or- 
der accuracy for a specific choice of slopes, without using 
PPM. 

It is true that third-order convergence is rarely seen in 
practice. In the context of Computational Fluid Dynam- 
ics (CFD), this is due to the arising of physical solutions 
(containing shocks or other discontinuities) which are 
just piecewise smooth. These discontinuities can prop- 
agate across the computational domain and the conver- 
gence rate is downgraded as a result in the regions swept 
away by the discontinuity front. A similar situation is 
encountered in black hole evolutions. The use of singu- 
larity avoidant slicing conditions produces a collapse in 
the lapse function. As it can be seen in Fig. [2j a steep 
gradient surface is formed (the collapse front) that prop- 
agates out as the grid points keep falling into the black 
hole. We will see that this results into a downgrade of 
accuracy in the regions close to the collapse front. 

Stability problems can also arise from the lack of res- 
olution of the collapse front, which is typically located 
around the apparent horizon. The reconstruction proce- 
dure can lead there to spurious oscillations, which intro- 
duce high-frequency noise in the simulation. In paper I, 
this problem was dealt with the use of standard slope lim- 
iters, which were crucial for the algorithm stability. In 
the present paper, although slope limiters are also dis- 
cussed for completeness, their use is not even required 
in any of the presented simulations. The new algorithm 
gets rid by itself of the high-frequency noise, even for the 
steep (but smooth) profiles appearing around the black- 
hole horizon. 

With all these simplifications, the proposed centered 
FV method can be interpreted just as an 'adaptive vis- 
cosity' generalization of the finite difference (FD) algo- 
rithms discussed before. Moreover, in the FV context, 
boundary conditions can be imposed in a simple way by 
the 'ghost point' technique. This allows one to avoid the 
complications related to the corners and edges treatment 
that usually appear in the FD context. 

The paper is organized as follows: we present in Sec- 
tion II a brief summary of the simplest FV methods. In 
Section III, the flux-splitting variant is considered, and 
we show how third-order space accuracy can be obtained 
by using just linear reconstruction. The resulting method 
is then tested for the one-dimensional (ID) black-hole in 
Section IV. Long term (up to 1000m) simulations are 
performed with a single numerical grid of a limited reso- 
lution, showing the efficiency of the algorithm. A conver- 
gence test is also performed, which confirms the predicted 
third-order accuracy in the outside region. The three- 
dimensional (3D) black-hole case is considered in Section 
V. A low resolution simulation is presented, showing the 
key role of controlling the trace of the extrinsic curvature 
in order to avoid numerical instabilities. This explains 



the advantage of using tr K as a primitive variable, like 
in the Conformal ADM (CADM) formalism Q. This 
explains also why a conformal decomposition was also 
required for obtaining robust 3D simulations in paper I, 
even when using FV methods [6|. 

For the sake of clarity, the more technical points: sta- 
bility analysis, time evolution algorithms and the full ex- 
plicit form of the equations, are described in Appendices 
A, B and C, respectively. 

II. CENTERED FINITE VOLUME METHODS: 
FLUX FORMULAE 

Let us consider the well known 3+1 decomposition of 
Einstein's field equations. The extrinsic curvature Kij 
is considered as an independent dynamical field, so that 
the evolution system is of first order in time but second 
order in space. Let us transform it into a fully first order 
system by considering also the first space derivatives of 
the metric as independent quantities. This requires ad- 
ditional evolution equations for these space derivatives, 
that can be obtained in the standard way by permuting 
space and time derivatives of the metric, that is 

d t (d k g a b) = d k (d t g a b) , (1) 

so that the resulting first order system will describe the 
same dynamics than the original second order one. 

In this first order form, Einstein's field equations can 
always be expressed as a system of balance laws Q . The 
evolution system can be written in the form 

9 t u + ^F fe (u) = S(u), (2) 

where both the Flux terms F and the Source terms S 
depend algebraically on the array of dynamical fields u, 
which contains the metric and all its first derivatives. 
The terms 'Fluxes' and 'Sources' come from the hydro- 
dynamical analogous of the system . 

The balance law form is well suited for FV discretiza- 
tion methods. The idea is to evolve the average of the 
dynamical fields u on some elementary cells, instead of 
evolving just point values like in the FD approach. The 
space discretization can be obtained by averaging ^ over 
an elementary cell and applying the divergence theorem 
to get: 

dt u + j> F fc dS k = S , (3) 

where the overlines stand for space averages. The eval- 
uation of partial space derivatives has been replaced in 
this way by that of surface integrals of the flux terms. 

Let us consider for simplicity the one-dimensional 
case. We can start from a regular finite difference grid. 
The elementary cell can then be chosen as the interval 
{ x i-i/2 i x i+i/i)i centered on the generic grid point Xi. 
The dynamical fields u can be modelled as piecewise 
linear functions in every cell (linear reconstruction, see 
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FIG. 1: Piecewise linear reconstruction of a given function. 
Numerical discontinuities appear at every cell interface (dot- 
ted lines) between the left and right values (arrows and dots, 
respectively). Note that the original function was monoton- 
ically decreasing: all the slopes are negative. However, both 
the left interface values (at i+3/2) and the right interface ones 
(at i — 3/2) show local extremes that break the monotonicity 
of the original function. 

Fig. [J), so that the average values Ui coincide with the 
point values u^. The corresponding (first-order accurate) 
FV discretization of is then given by 

< +1 = < - ^ [ F* +1/2 - FU/2 ] + Af Si . (4) 

We will restrict ourselves to these linear reconstruction 
methods in what follows. 



Flux formulae 

The generic algorithm Q requires some prescription 
for the interface fluxes Py-1/2 ■ A straightforward calcu- 
lation shows that the simple average 

Fi+i/2 = \ (Fi + F i+1 ) (5) 

makes ((H) fully equivalent to the standard second order 
FD approach. As it is well known, this choice is prone to 
developing high-frequency noise in presence of steep gra- 
dients, like the ones appearing in black hole simulations. 
For this reason, artificial viscosity terms are usually re- 
quired in order to suppress the spurious high-frequency 
modes [J]. 

We will consider here more general flux formulae, 
namely 

F i+ i/2 = f(u L , u R ), (6) 

where Ul, ur stand for the left and right predictions for 
the dynamical field u at the chosen interface (arrows and 
dots, respectively, in Fig. [1]). In the (piecewise) linear 
case, they are given by 

u L = Ui + 1/2 o i Ax u R = u l+1 -1/2 a i+1 Ax , (7) 



where cr, stands for the slope of the chosen field in the 
corresponding cell. 

A sophisticated choice is provided by the 'shock- 
capturing' methods (see Ref. Q for a review). The idea 
is to consider the jump at the interface as a physical one 
(not just a numerical artifact). The characteristic de- 
composition of (the principal part of) the system is then 
used in order to compute some physically sound inter- 
face Flux. These advanced methods have been common 
practice in Computational Fluid Dynamics since decades. 
They were adapted to the Numerical Relativity context 
fifteen years ago 0] , for dealing with the spherically sym- 
metric (ID) black-hole case. They are still currently used 
in Relativistic Hydrodynamics codes, but their use in 3D 
black hole simulations has been limited by the compu- 
tational cost of performing the characteristic decomposi- 
tion of the evolution system at every single interface. 

More recently, much simpler alternatives have been 
proposed, which require just the knowledge of the charac- 
teristic speeds, not the full characteristic decomposition. 
Some of them have yet been implemented in Relativistic 
Hydrodynamics codes Maybe the simplest choice is 
the local Lax-Friedrichs (LLF) flux formula 

f(u L , u r ) = ~[F l +Fr + c(u l - Ur) } , (8) 

where the coefficient c depends on the values of the char- 
acteristic speeds at the interface, namely: 

c = max( X L , X R ) , (9) 

where A is the spectral radius (the absolute value of the 
biggest characteristic speed). 

When comparing the LLF choice ([5]) with the centered 
FD one (O, we can see that the supplementary terms 
play the role of a numerical dissipation. In this sense, a 
much more dissipative choice would be 



which corresponds to (a piecewise linear generalization 
of) the original Lax-Friedrichs algorithm. Note that in 
any case the values of the dissipation coefficients are pre- 
scribed by the numerical algorithms: they are no arbi- 
trary parameters, like in the FD case. 



III. FLUX SPLITTING APPROACH 

In the flux formulae approach ([6]), the information 
coming from both sides is processed at every interface, 
where different components are selected from either side 
in order to build up the flux there. We will consider 
here an alternative approach, in which the information 
is processed instead at the grid nodes, by selecting there 
the components of the flux that will propagate in either 
direction (flux splitting approach) [5(. 
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The flux-splitting analogous of the original LLF for- 
mula ([HI can be obtained by splitting the flux into two 
simple components 

F l ± = F i ± A, Ul , (11) 

where A will be again the spectral radius at the given 
grid point. Each component is then reconstructed sep- 
arately, leading to one-sided predictions at the neighbor 
interfaces. The final interface flux will be computed then 
simply as 

Fi+i/2 = \ {F+ + F-). (12) 

This method can also be expressed as a modified LLF 
flux formula, namely 

f(uL , u R ) = - [ F L + F R + X L u L - X R u R ] . (13) 

The main difference between the original LLF flux for- 
mula © and the flux-splitting variant (fT3")) is that in the 
last case there is a clear-cut separation between the con- 
tributions coming from either the left or the right side 
of the interface, as it can clearly be seen in IT2"]) . In this 
way, one has a clear vision of the information flux in the 
numerical algorithm. The information from F + compo- 
nents propagates in the forward direction, whereas the 
one from F~ components propagates backwards. This 
simple splitting provides in this way some insight that 
can be useful for setting up suitable boundary conditions. 
Moreover, it opens the door to using different slopes for 
the reconstruction of each flux component. We will see 
below how to take advantage of this fact in order to im- 
prove space accuracy. 

Third order accuracy 

As it is well known, the use of a consistent piecewise- 
linear reconstruction results generically into a second- 
order space accuracy. A convenient choice is given by 
the centered slope 

° C = 2^ ( Ul+1 ~ Ui ~^' ^ 

This is a good default choice (Fromm choice @), leading 
to reliable second-order accurate algorithms . 

More general second-order algorithms can be obtained 
by replacing the centered slope a by any convex average 
of the left and right slopes, 

cr L = (m - Ui-i)/Ax , a R = (uj+i - Ui)/Ax . (15) 

In some applications, however, second order accuracy is 
not enough. The leading (third order) error is of the dis- 
persion type, affecting the numerical propagation speeds. 
In the FD approach, this can be improved by using a 
fourth-order-accurate algorithm in combination with a 



fourth-order artificial dissipation term (which constitutes 
itself the leading error term). The resulting combination 
is third-order accurate. 

In the standard FV approach, the standard way of get- 
ting (piecewise) third-order accuracy would be instead to 
replace the piecewise linear reconstruction by a piecewise 
parabolic one. The prototypical example is provided by 
the well known piecewise parabolic methods (PPM) . The 
main complication of this strategy is that node values 
would no longer represent the cell averages of a given dy- 
namical field. This would increase the complexity of the 
reconstruction process and the computational cost of the 
resulting algorithm. 

There is a much simpler alternative, which takes ad- 
vantage of the Flux splitting (fTTj) . The idea is to consider 
the resulting one-sided components F^ as independent 
dynamical fields, each one with its own slope. The sur- 
prising result is that the choice 

* + = l* L + 2 -a R , a- = lo* + \oX (16) 

leads, after the recombination (fl2|) . to a third-order accu- 
rate algorithm. The coefficients in (| are unique: any 
other combination leads just to second-order accuracy. 

Note that we are getting in this way third-order ac- 
curacy with a piecewise linear reconstruction (see the 
convergence test in Fig. [5] for a confirmation). This 
important result seems to be a peculiarity of the Flux- 
splitting approach. In order to better understand it, let 
us suppress for a moment the lambda terms in (|lHll3p . 
A straightforward calculation shows that, when using 
the slopes (|16l) , the resulting algorithm coincides exactly 
with the standard fourth-order-accurate FD algorithm. 
Adding the lambda terms improves the stability of the 
algorithm at the price of downgrading the space accuracy 
to third order. This is precisely the same effect that the 
Kreiss-Oliger dissipation terms produce in the FD case. 
This confirms our result and suggests the interpretation 
of the algorithm (|1HI13[) as providing an adaptive gener- 
alization of the standard dissipation terms. 

IV. THE ID BLACK HOLE 

As a first test, let us consider the Schwarzschild Black 
Hole in spherical coordinates. We will write the line ele- 
ment in the 'wormhole' form: 

ds 2 = -( tanhrj f dt 2 + 4m 2 ( coshij/2 ) 4 ( drf + dtt 2 ) , 

which can be obtained from the isotropic form by the 
following coordinate transformation 

r = to/2 exp ( r] ) . (18) 

The wormhole form (| L7[) exploits the presence of a min- 
imal surface (throat) at n — 0. It is manifestly invariant 
by the reflection isometry 

r\ «-> — T] , (19) 
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FIG. 2: Long-term FV simulation of a ID black hole, with 
a single mesh of 120 gridpoints. The evolution of the lapse 
is shown up to 1000m, in intervals of 50m (solid lines). The 
dotted lines correspond to lm, 3m, 5m and 25m. Note that 
the plots tend to cumulate at the end, due to the exponential 
character of the grid, as given by (|18[) . No slope limiters have 
been used in this simulation. 



so that the numerical simulations can be restricted to 
positive values of rj. The isometry (|19[) provides a very 
convenient boundary condition at the throat. Moreover 
(HHJ) implies 



dr = r dr] 



(20) 



so that an evenly spaced grid in r\ corresponds to a ge- 
ometrically increasing spacing in r. We can perform in 
this way long term simulations with a single grid of a 
limited size, as we will see below. This also allows to 
apply the standard boundary conditions in FV methods: 
two 'ghost' points are added by just copying the near- 
est neighbor values (or their time variation) for every 
dynamical field. The separation between incoming and 
outgoing information is automatically performed by the 
flux-splitting algorithm, so that boundary points are not 
special in this respect. 

The simulations are performed with a spherically sym- 
metric version of the Z3 formalism as detailed in 
Appendix C. The free parameter n, governing the cou- 
pling with the energy constraint, is taken with unit value 
by default, but other similar values can be taken without 
affecting significatively the results, like n = 4/3, which 
corresponds to the CADM case [15(. Regarding gauge 
conditions, we are using the generalized harmonic pre- 
scription for the lapse [l6| 



{d t - Cp) a = -f a 2 trK 



(21) 



with zero shift (normal coordinates). We take a constant 
(unit) value of the lapse as initial data. We can see in 
Fig. [2] the evolution of the lapse in a long-term simulation 
(up to 1000m). We have chosen in this case / = 2/ a 
(corresponding to the 1+log slicing), but similar results 



FIG. 3: The evolution of the propagation speed is shown up 
to 100m, in intervals of 10m, for the same simulation as in 
Fig. The maximum values are clearly seen to decrease in 
time. Note the exponentially decreasing tail, as a result of 
the choice of the radial coordinate. 



can be obtained with many other combinations of the 
form 



f = a + b/c 



(22) 



where a and b are constant parameters. 

Note that no slope limiters have been used in the sim- 
ulation shown in Fig. O This can seem surprising at the 
first sight, but it can be better understood by looking 
at the propagation speed profiles shown in Fig. [3] The 
maximum propagation speed values decrease with time, 
due to the lapse collapse in the black hole interior re- 
gion. This happens because the initial speed profile is 
exponentially decreasing with the chosen radial coordi- 
nate. The same decreasing arises for gauge speed. As a 
result, the Courant stability condition becomes less and 
less restrictive as the simulation proceeds, allowing us to 
take bigger timesteps. We have preferred instead to keep 
the initial timestep for the sake of accuracy. As far as all 
derivative terms get multiplied by At in the algorithm 
(H|), this gives us an extra safety factor that allows us to 
avoid using slope limiters. 

As an accuracy check, we monitor the mass func- 
tion |17| , which is to be constant in space and time for the 



Schwarzschild case, independently of the coordinate sys- 
tem. In Fig. 31 we compare (the Li norm of) the errors 
in the mass function between a third-order FV simula- 
tion (without slope limiters) and the corresponding FD 
simulation (including a fourth order dissipation term like 
the one in ref. [3] with e = 0.015). We see that the FD 
method shows bigger errors at late times. One can argue 
that the leading error in the FD simulation is given by the 
dissipation terms, so that one can modify the result by 
lowering the numerical dissipation coefficient. However, 
lowering the viscosity coefficient used in Fig. [H would re- 
sult into a premature code crashing, like the one shown 
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FIG. 4: Time evolution of the error in the mass function 
(logarithm of the L2 norm) for three different numerical al- 
gorithms. The strictly fourth-order FD method, without ex- 
tra dissipation terms, is the most accurate as expected, but 
crashes after a short time (measured in units of m). The 
other two algorithms (third-order accurate) get similar errors 
at early times, but the FV one performs much better in the 
long term than the FD with standard Kreiss-Oliger dissipa- 
tion. The dissipation coefficient has been taken as low as 
allowed by code stability (see the text). All simulations were 
obtained with a single mesh of 120 gridpoints and using the 
1+log slicing prescription. 



in the Figure for a strictly fourth order FD run, without 
the artificial dissipation term. 

We can understand the need for dissipation by looking 
at the sharp collapse front in Fig. [2] We know that this 
is not a shock: it could be perfectly resolved by increas- 
ing the grid resolution as needed. In this way we can 
actually get long-term ID black hole simulations, with 
a lifetime depending on the allowed resolution. This 
'brute force' approach, however, can not be translated 
into the 3D case, where a more efficient management of 
the computational resources is required. This is where 
dissipation comes into play, either the numerical dissi- 
pation built in FV methods or the artificial one which 
is routinely added to fourth-order FD methods. Dissi- 
pation is very efficient in damping sharp features, corre- 
sponding to high-frequency Fourier modes. As a result, 
the collapse front gets smoothed out and can be resolved 
without allocating too many grid points. However, the 
more dissipation the more error. In this sense, Fig. [4] 
shows that adaptive viscosity built in the proposed FV 
method provides a good compromise between accuracy 
and computational efficiency. 

Note that the error comparison is independent of the 
selected resolution. This is because the two stable meth- 
ods in Fig. 0] are of third order accuracy, as confirmed by 
the local convergence test shown in Fig. [5] (solid line, cor- 
responding to t — 10 to). In the long term, however, large 
errors develop around the collapse front, downgrading the 
local convergence rate in the neighbor regions (dashed 



FIG. 5: Local convergence evolution for the mass function in 
a ID black hole simulation. We can see the predicted third- 
order accuracy, when using the proposed slopes (|16p . around 
t = 10m (solid line). At t = 100m (dashed line), we yet 
see the downgrade in the regions around the collapse front 
(the apparent horizon position is marked with a circle). As 
the collapse front propagates (dotted line, corresponding to 
t — 400 m), we can see the growth of the affected regions, 
specially the one behind the front. 



and dotted lines in Fig. corresponding to t = 100 m 
and t = 400m, respectively). This can not be seen as 
a failure of the algorithm properties, but rather as con- 
sequence of large errors in a highly non- linear context. 
This also shows that in simulations oriented to compute 
gravitational wave patterns (not the case of this paper), 
the waveform extraction zone must be safely located, 
away both from the outer boundary and from the col- 
lapse front. 



V. PRELIMINARY 3D RESULTS 

The ID algorithm ([4]) can be easily adapted to the full 
three-dimensional (3D) case: 



U {ijk} ~ U {tjk} 



I F {i+l/2ife} ^ F {i-l/2jk} 1 

At r -pvy _ Tpjf 

Ay 1 {ij + l/2fc} F {ij-l/2k} 

At r _. 



Az 
At S{ijk} 



{ij fe+l/2} r {ij fc-1/2} 



(23) 



The structure of (|23| suggests dealing with the 3D prob- 
lem as a simple superposition of ID problems along every 
single space direction. The stability analysis in Appendix 
A can then be extended in a straightforward way, show- 
ing that the strong stability requirement leads to a more 
restrictive upper bound on the timestep (in our case, us- 
ing a cubic grid, this amounts to an extra 1/3 factor). 



7 




-0.4 1 1 1 1 1 1 1 1 

1 2 3 4 5 6 7 

FIG. 6: Plot of the trace of the extrinsic curvature at t = 12m 
for a low resolution simulation. The dotted line corresponds 
to the trace obtained by contraction from the individual com- 
ponents Kij . The solid line is the same quantity computed di- 
rectly as a primitive variable. The big difference corresponds 
to the transition between the collapsed and uncollapsed re- 
gions, where the lapse shows a steep profile 

In cartesian-like coordinates, it is not so easy to take 
advantage of the reflection isometry (| 10[) . For this rea- 
son, we will evolve both the black-hole exterior and the 
interior domains. We can not use the r\ coordinate for 
this purpose, because the symmetry center would corre- 
spond to r\ — > oo. We will take instead the initial space 
metric in isotropic coordinates, namely 

dl 2 = (1 + — ) 4 Sa dx'dx 3 . (24) 
2r 

We will replace then the vacuum black-hole interior by 
some singularity-free matter solution. To be more spe- 
cific, we will allow the initial mass to have a radial de- 
pendence: m = m(r) in the interior region. This allows 
to match a scalar field interior metric to ('stuffed 
black-hole' approach [HI). The price to pay for using a 
regular metric inside the horizon is to evolve the matter 
content during the simulation: we have chosen the scalar 
field just for simplicity. 

We have performed then a low-resolution simulation 
(Ax = 0.1m) in order to monitor the errors in tr K, which 
determines the evolution of the lapse. We see in Fig. [5] 
the comparison between the trace computed by contract- 
ing the individual components (dotted line) and an 
auxiliary variable K which is evolved by using the an- 
alytical equation for tr K (solid line). The difference is 
striking, even at the early time of the plot (t = 12m). 
Note the negative peak in the computed tr K , which will 
produce a spike in the lapse leading to a premature code 
crashing. 

This behavior could be somehow anticipated from our 
previous ID simulations. The plots shown in Fig. [5] cor- 
respond to the mixed indices equations displayed in Ap- 
pendix C. We have performed for comparison the same 



simulations with 'downstairs' indices and the results look 
different. We actually double-checked both codes before 
realizing that just raising one index can make a difference 
at a given resolution. Of course, in ID we can always in- 
crease resolution at will and verify that the two results 
get close enough. But this would be prohibitive in 3D, 
at least for single-grid simulations. Moreover, in 3D we 
have the additional difficulty of modelling curved features 
in a Cartesian grid. In the spherical case, the worst sit- 
uation shows up along the main diagonal, precisely the 
view shown in Fig. 

These considerations can explain why the CADM for- 
malism, which actually uses iv K as a primitive variable, 
has shown to be more robust even in single-grid simula- 
tions. This also explains why the use of a conformal de- 
composition was crucial in the 3D simulations performed 
with the old (non-covariant) Bona-Masso formalism in 
paper I, which used shock-capturing methods. The Z3 
formalism can be interpreted as a covariant version of 
the same, but our results strongly suggest that the key 
element for robustness is not covariance but the use of a 
conformal decomposition. 

As a final remark, let us focus on the boundary con- 
ditions implementation. The 3D FV algorithm (|23|) al- 
lows to apply the ghost point technique exactly in the 
same way as in the ID case: by just copying (the time 
variation of) all the quantities from the neighbor inte- 
rior point. There is no need for any special treatment 
for corners or vertices. Moreover, the simple FV meth- 
ods presented here do not require the explicit use of the 
characteristic decomposition, not even at the boundaries. 
In spite of these simplifications, the robust stability test 
for the combined initial-boundary problem gives results 
equivalent to the ones obtained with maximally dissipa- 
tive boundary conditions in a finite difference context 
(see Appendix B in Ref. for details). 
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Appendix A: Stability and Monotonicity 

Let us assume that (the principal part of) the evo- 
lution system is strongly hyperbolic. This means that, 
for any chosen direction, we can express the system as 
a set of simple advection equations for the characteris- 
tic variables (eigenfields). In order to verify the stability 
properties of the proposed algorithms, it will be enough 
to consider a single advection equation with a generic 
speed v. The corresponding Flux will be given then by 

F(u) = v u . (25) 
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We will consider in the first place the first-order accu- 
rate approximation, obtained by a piecewise constant re- 
construction (zero slope). The corresponding discretiza- 
tion can be obtained by replacing the prescription Q12p 
into the general expression ([4]). The result is the linear 
three-point algorithm: 



,n+i 



~Ax [ 2 (Al+1 
1 



Vi+l) <+l 



+ - (Aw+M«i-i-Ai<] -(26) 

Allowing for the fact that A is chosen at every point as 
the absolute value of the maximum speed, we can see 
that all the u n coefficients are positive provided that the 
Courant stability condition 



, At 
A — < 1 
Ax ~ 



(27) 



is satisfied. Note however that a more restrictive condi- 
tion is obtained in the three-dimensional case, where we 
must add up in (|26p the contributions from every space 
direction. 

As it is well known, the positivity of all the coefficients 
ensures that the algorithm is monotonicity-preserving, 
so that spurious numerical oscillations can not appear. 
This implies stability, but the converse is not true, as it 
is well known. Let us remember at this point that the 
centered FD discretization could be recovered from 
simply by setting A to zero, although we would lose the 
monotonicity property in this way. 

The monotonicity properties of the piecewise constant 
reconstruction are not ensured in the piecewise linear 
case. We can clearly see in Fig. [1] that monotonicity 
problems can arise in steep gradient regions. The reason 
is that either the series of left {u L } or right {u R } in- 
terface predictions can show spurious peaks which where 
not present in the original function. In the case of the 
centered slope ([5]), a detailed analysis shows that this 
will happen at a given interface only if the left and right 
slopes differ by a factor of three or more. This gives a 
more precise sense to the 'steep gradient' notion in the 
centered slopes case. 

The natural way to remedy this is to enforce that both 
(left and right) interface predictions are in the interval 
limited by the corresponding left and right point values 
(interwinding requirement). This amounts to using the 
'limited' slopes 



lira 



minmod( 2a 1 



2a l 



(28) 



where a is the default slope at the given cell. This inter- 
winding requirement is not enough, however, to ensure 
the positivity of all the coefficients in the resulting algo- 
rithm. A detailed analysis shows that an extra factor in 
the Courant condition would be required for monotonic- 
ity in this case: 



A^<l/2 
Ax ~ 1 



Note however that we are analyzing here the elementary 
step (|4|. This is just the building block of the time evo- 
lution algorithm. The exact stability and monotonicity 
limits for the time step would depend on the specific 
choice of the full time evolution algorithm [4] , which will 
be described in Appendix B. 

A word of caution must be given at this point. It is well 
known that the monotonicity results hold only for strictly 
Flux-conservative algorithms. This is not our case: the 
Source terms play an important physical role. Of course, 
these terms do not belong to the principal part, so that 
positivity of the Flux terms ensures some strong form 
of stability. Nevertheless, one must be very careful with 
the physical interpretation, because the first-order con- 
straints |T]) preclude any clear-cut isolation of the Source 
terms. This makes the analogy with Fluid Dynamics only 
approximative and the use of the slope limiters a risky 
matter: we could be removing in the Flux part some fea- 
tures that are required to compensate something in the 
Source part. Our experience is that, at least for smooth 
profiles, more robust numerical simulations are obtained 
when the slope limiters are switched off. The high fre- 
quency modes are kept under control by the numerical 
dissipation built in the proposed FV methods. 



Appendix B: Time accuracy 

The simple step Q is only first-order accurate in time, 
and this fact is not changed by any of the space accuracy 
improvements we have considered up to now. The stan- 
dard way of improving time accuracy is by the method 
of lines (MoL, see refs. [12| |4j). The idea is to consider 
(HI as a basic evolution step 



E( u" , At 



(30) 



in order to build higher order algorithms. A convenient 
choice for these time evolution algorithms is provided the 
standard Runge-Kutta methods [l3| (see also Q). For 
instance, second order accuracy can be obtained in two 
steps: 

u* = E{u n , At ) u n+l = - [u n +E(u*, At)}, (31) 

and third-order time accuracy with one more intermedi- 
ate step: 

u** = - A u n + - A E{ u* , At) 



~ 3 1 3 



«" ' = - u- + - E( u** , At ) 



(32) 



(29) 



Note that the positivity of all the coefficients in 
ensures that the monotonicity property of the basic step 
([30)1 will be preserved by the resulting strong-stability- 
preserving (SSP) algorithm. This interesting property 
comes at the price of keeping the upper limit on At that 
is required for the monotonicity of the basic step. This 
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is a clear disadvantage with respect to the case in which 
the standard FD approach is being used for space dis- 
cretization, in which one is only limited by plain sta- 
bility, not monotonicity. Then, there are Runge-Kutta 
algorithms (with non-positive coefficients) that alow to 
take At larger than the one required by the standard 
Courant condition 

Conversely, second order Runge-Kutta algorithms like 
(I3ip are unstable when used in combination with FD 
space discretization, unless artificial dissipation is added 
in order to recover stability (not just monotonicity) 
This is why FD simulations currently use at least a third- 
order time evolution algorithm. 

Appendix C: Z3 evolution equations 

The Z3 evolution system [l4|, HB| is given by: 

(d t - Cp) jij = -2aKij (33) 
{d t - C fj )K l3 =-V i a j + a[ R tJ + V t Z + VjZi 

- 2 Kfj + trK Kij - Sij + ^ ( trS + (n - 1) r ) 7y ] 

- ^ a[trR + 2 X7 k Z k 

+4 tr 2 K - tr(K 2 ) - 2 Z k a k /a] 7^ (34) 
(d t - C )Zi=a [Vj (K* - 5i j trK) - 2K i j Z j - St] , 

(35) 

where n is an arbitrary parameter governing the coupling 
of the energy constraint. 

The fully first-order version can be obtained in the 
standard way, by introducing the additional fields 

D k ij = - d k Jij ■ (36) 
Note that the ordering constraint reads 

d r Dkij = dk D ri j , (37) 

which is no longer an identity for the first order system. 
As a consequence of this ordering ambiguity of second 
derivatives, the Ricci tensor term in (the first order ver- 
sion of) the evolution equation (|34j) can be written in 
many different ways. Then, an ordering parameter Q can 
be introduced (l5| . so that the parameter choice £ = +1 
corresponds to the standard Ricci decomposition 

^ -^i? = ij T kj -\- r r fcY ij r r {Y kj (38) 

whereas the opposite choice £ = — 1 corresponds instead 
to the decomposition 

Rij = — &k D k ij + 9(i Yj) k k — 2D r rk Dkij 

_i_ A firs j~) _ p p rs _ p prfc CiQ\ 



which is most commonly used in Numerical Relativity 
codes. We can then consider the generic case as a linear 
combination of (|38|) and (|39| . 

In the spherically symmetric vacuum case, the first or- 
der version of the system (I33ti34[) is free of any ordering 
ambiguity. It can be written as 

dtirr = -2aj rr K r r , d t jee = -2ajeeK e (40) 

d t K r r + d r [a Y r (A r + (2 - n) D 9 ° - (2 - n/2) Z r )] = 
a [{K/f + (2 - n) K r r K e e - (n/2) [K e e f 
-Y r D r r (At + (2 - n) Dg" + {n/2 - 2) Z r ) 
+7 rr D 9 9 {{2-n)A r -{2- in/2) Dg 8 - n Z r ) 
~Y r (2~n)A r Z r -(n/2)-f ee ] (41) 

8tK e ° + d r [a 7 rr {(l-n)D e e + {n/2) Z r )] = 
a[(l-n)K/K d s +{2 -n/2) (K g e ) 2 
-Y r D r r ((l-n)D g e + (n/2)Z r ) 
+ 7 rr Dg 9 ((2-n)Z r -(2- 3n/2) D g ") 
- ni rr A r {D e ° - Z r ) + (l-n/2) 1 ee ] (42) 
d t Z r + d r [2aK e 6 ] = 

2 a [D e e (K r r - K 9 e ) + A r K e e - K r r Z r ] (43) 

d t D r r + d r [a K T r \ = 0, d t Dg 6 + d r [a K g e ] = 0, (44) 

where we are using normal coordinates (zero shift). The 
slicing condition (f2Tj) can be written as 

d t a = -a 2 f trK , d t A r + d r [a f trK] = . (45) 



The mass function can be defined for spherically sym- 
metric spacetimes as [l7j 

2M = Y [ l-g ab d a Yd b Y } , (46) 



where Y stands for the area radius. In spherical coordi- 
nates we get 

2M(t, r) = { 1 + lee [(Kg 8 ) 2 - ^(Dg 6 ) 2 } } . (47) 



The mass function has a clear physical interpretation: 
it provides the mass inside a sphere of radius r at the 
given time t. It follows that M(t, r) must be constant 
for the Schwarzschild spacetime, no matter which coor- 
dinates are being used. This provides a convenient accu- 
racy check for numerical simulations. 
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